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Abstract Defects are believed to play a fundamental role in the supersolid 
state of 4 He. We report on studies by exact Quantum Monte Carlo (QMC) 
O ■ simulations at zero temperature of the properties of solid 4 He in presence of 

a many vacancies, up to 30 in two dimensions (2D). In all studied cases the 

crystalline order is stable at least as long as the concentration of vacancies 
is below 2.5%. In the 2D system for a small number, n v , of vacancies such 
defects can be identified in the crystalline lattice and are strongly correlated 
with an attractive interaction. On the contrary when n v > 10 vacancies in 
the relaxed system disappear and in their place one finds dislocations and 
a revival of the Bose-Einstein condensation. Thus, should zero-point motion 
defects be present in solid 4 He, such defects would be dislocations and not 
vacancies, at least in 2D. In order to avoid using periodic boundary conditions 
we have studied the exact ground state of solid 4 He confined in a circular 
region by an external potential. We find that defects tend to be localized in 
an interfacial region of width of about 15 A. Our computation allows to put 
as upper bound limit to zero-point defects the concentration 3 x 10" 
2D system close to melting density. 
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1 Introduction 

Supersolidity is an intriguing state of matter in which spatial order, typical 
of the solid phase, and off-diagonal long range order [T] (ODLRO), which 
characterizes the Bose-Einstein condensation (BEC) phenomena, are simul- 
taneously present, implying some form of superfluid properties even for the 
solid phase. This striking idea of a superfluid solid was proposed long ago 
[3] and solid 4 He was early recognized as the natural candidate to display such 
a counterintuitive coexistence of orders. This topic attracted the attention 
of many physicists 4 before going out of attention in the absence of exper- 
imental evidence. As a consequence of the recent discovery of non classical 
rotational inertia (NCRI) [5] in solid 4 He, one of the expected manifestations 
of supersolidity [5J , the possible existence of a supersolid phase has recently 
gained once more the attention of the scientific community, making solid 4 Hc 
systems the subject of many experimental and theoretical studies (See Refs. 
[71151 151110] for recent reviews). 

There is strong evidence that defects play an important role in NCRI. 
In fact, experiments show that usually NCRI is strengthened by increasing 
disorder in the crystal, and microscopic simulation studies agree on the fact 
that an ideal perfect crystal does not show ODLRO JTTIITS] even at T = 
K |13l l9"] . By ideal perfect crystal we mean a crystalline solid extended to all 
space with one atom per lattice site. A key question is if bulk solid 4 He is an 
ideal perfect crystal or not. For instance, if zero-point independent vacancies 
were present, the crystal would still be perfect because vacancies are delo- 
calized defects, but not ideal because there would be less than one atom per 
lattice site (this is also called an incommensurate solid) . Notice that the ear- 
lier proposal [3]|3] for the supersolid state assumed the presence of zero-point 
defects in the ground state of solid 4 He in the form of ground state vacancies. 
Presently there is a prevalent view that supersolidity is an extrinsic property 
of solid 4 He. However not all share this view jTHf5|IT5] and present experi- 
ments cannot rule out the presence of intrinsic defects [To¥l~7] . Moreover, it 
is still not clear which kind of disorder contributes to, or is responsible for, 
the anomalous properties observed in solid 4 He. Many different defects have 
been considered, ranging from vacancies to grain boundaries, passing from 
dislocations to quantum glasses; but none of the proposed models seems able 
to capture the whole phenomenology of supersolidity in solid 4 He [TUJ. Thus, 
independent of the very relevant issue of intrinsic or extrinsic defects, the 
knowledge of the properties of defects in this system is one of the main goals 
of the present theoretical investigations on supersolidity in this system. 

From the theoretical point of view the method of choice to investigate a 
strongly interacting quantum system like solid 4 He is Quantum Monte Carlo 
(QMC) simulations. Nowadays a number of studies on defected solid 4 He sys- 
tems studied by means of QMC techniques are present in literature [1811191 
|^|^|^|^[^I^I^I^I^I3TTO^ . Crystals with de- 

fects can be simulated with an appropriate choice of the simulation box (SB) 
and the number of particles such that, combined with the use of periodic 
boundary conditions (PBC), the MC sampling is constrained to configura- 
tions which host the desired defect. When PBC are not able to stabilize the 



defect, one can constrain (or fix) the degrees of freedom of a number of atoms 
surrounding the defect of interest [30]. These constraints on the configura- 
tional space, if judiciously implemented, usually do not prevent the study of 
the physical properties of the (defected) system; as an example, the binding 
energy of a 3 He atom to dislocation cores turned out to be in agreement 
with experimental data |31) . A distinct issue is the investigation of the na- 
ture (defected or not) of the ground state of the extended system; in this 
case, such strategies are not straightforwardly convincing and some doubt 
remains on the level of confidence to give to the results: beside the finite 
size effects, one also has to check to what extent the results are influenced 
by the choice of boundary conditions. In fact, the constraint imposed by the 
PBC on the configurational sampling of a crystalline solid is both stronger 
and subtler than what is usually believed, due to commensurability effects 
between crystal unit cell and SB |40) . 

In the present paper we address two topics. The first is to study if it is 
true that vacancies are not a viable route to supersolidity because multiple 
vacancies coalesce and lead to phase separation [33] ■ The second topic is to 
study what we learn from quantum simulations by studying large systems 
without using PBC. On the first topic we review and extend a systematic 
study of many vacancies (up to 30) in two-dimensional (2D) solid 4 He based 
on exact T — K quantum simulations. A pure 2D solid is of interest as 
a simple model for out of registry solid 4 He adsorbed on planar substrates, 
and new experimental investigations |39j are under way to answer the ques- 
tion of whether a supersolid state is present also in adsorbed 4 He. In no 
case do we find that the crystalline order is unstable, at least as long as 
the concentration of vacancies is below 2.5%. Multiple vacancies were found 
highly correlated with a tendency to form linear structures. Actually, when 
n v > 10 vacancies transform themselves into dislocations. This resembles the 
behavior of multiple vacancies in classical solids at low temperature |41j . In 
our quantum crystal such dislocations are very mobile, allowing exchange of 
particles across the system; in fact, we find ODLRO in the one-body density 
matrix p\ computed for the finite simulated system. In presence of disorder, 
this gives the possibility of establishing a well defined phase [42 , at least 
locally, in the extended system. Should the exact ground state of 4 He in 2D 
not be an ideal perfect crystal, our results imply that the zero-point disor- 
der does not consist of vacancies but of dislocations and that the extended 
system could be supersolid. 

As second topic we study systems free of PBC by confining the 4 He atoms 
in a finite region of space by an external potential. The main aim is to inves- 
tigate the true nature of the ground state, i.e. to determine if defects restrict 
themselves in the intcrfacial region near to the confining potential, triggered 
by the mismatch between the confining geometry and the triangular lattice 
of 2D solid 4 He, or if defects also permeate the inner region. We give evidence 
that once the equilibrium is reached, the defects, even if initially placed in the 
center of the system, essentially restrict themselves in the interfacial region, 
leaving the inner region a regular triangular lattice. Moreover, we find that 
the inner crystal is compatible with the ideal crystal simulated in a periodi- 
cally repeated box (i.e. in standard simulations with PBC). Given the size of 



the studied systems, this gives an upper bound on the concentration of any 
ground state defects, which must be below 3 x 10~ 3 . 

The paper is organized as follows: Sec. [5] deals with the exact T = 0K 
shadow path integral ground state (SPIGS) method and the definition of 
quantities used in measuring the crystal disorder. Details on the simulations 
and our results on the periodically repeated crystal are presented in Sec. [3] 
Sec. @] contains a discussion on the issue of zero-point defects and on the 
influence of using PBC. Sec. \5\ contains details on the simulations and results 
on the confined 2D crystal. Conclusions are given in Sec. El 



2 Simulation details 

The SPIGS method: We study solid 4 He via exact simulation methods based 
on the Shadow Path Integral Ground State (SPIGS) [2Tj|22] method. SPIGS 
is an extension of the Path Integral Ground State (PIGS) [J3] method. The 
aim of PIGS is to improve a variationally optimized trial wave function tpT 
by constructing, in the Hilbert space of the system, a path which connects 
the given ipx with the exact lowest energy wave function of the system, ipo, 
constrained by the choice of the number of particle N, the geometry of the SB, 
the boundary conditions and the density p. During this "path" , the correct 
correlations among the particles arise through the "imaginary time evolution 
operator" e~ rH , where H is the Hamiltonian operator. For a large enough r, 
an accurate representation for the lowest energy state wave function is given 
by ip T = e~ TH ipT, which can be written analytically by discretizing the path 
in imaginary time; this maps the quantum system into a classical system of 
open polymers [43]. An appealing feature peculiar to the PIGS method is 
that, in ip T , the variational ansatz acts only as a starting point, while the 
full path in imaginary time is governed by e~ rH , which depends only on 
the Hamiltonian operator. We have recently shown that ip T does not need 
to be a variational optimized wave function: PIGS results for large enough r 
are unaffected by the choice of ipT both in the liquid and in the solid phase 
[T3"ll4"4] thus providing an unbiased exact T = K QMC method. Within 
SPIGS a shadow wave function (SWF) [4"51l4l>] is taken as ifo and this choice 
was shown to greatly accelerate convergence to ipo- Another feature of the 
SPIGS method is that it recovers the solid phase via a spontaneously broken 
translational symmetry [2"Tll2"2"] . like in Path Integral Monte Carlo (PIMC), 
and there is no constriction on the atomic positions, so it is particularly 
useful in studying crystals with defects. 

Measuring crystal order and fluctuations: We check the presence of crys- 
talline order by monitoring the static structure factor for the presence of 
Bragg peaks and, via a Delaunay triangulation [35] (DT) of the sampled 
configurations, we compute the particle coordination number in order to es- 
timate the amount of local disorder in the system [?7] . In an ideal perfect 2D 
triangular crystal, with the atoms in their equilibrium positions, each atom 
is linked to 6 other atoms in the DT. Atoms with coordination number not 
equal to 6 are then a measure of fluctuations and of local disorder in the 
crystal. In periodically repeated systems in 2D one proves a conservation law 



for the coordination numbers: 
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£(6-*)JVi = £;(*-6)iVi (1) 

i=3 i+7 

where Ni is the number of i— coordinated atoms, so that we can consider 
only Ni with i < 6. Since it is always verified that N% = 0, the quantity 
Xd = 2A^4 + N§ is usually taken as an estimate of local disorder in the 
system [48]. However, even in an ideal (defect- free) 2D quantum crystal the 
coordination is 6 only on the average: atoms are not always 6-fold coordinated 
due to their large zero-point motion. Then, a more useful quantity measuring 
the net amount of disorder in a crystal with n v vacancies, for instance, is 
the difference between the observed Xd and that of the corresponding ideal 
perfect crystal: 

X d = (2iV 4 + 7V 5 )„„ - (27V 4 + Ar 5 )n„=o. (2) 

We have studied also the orientational order parameter Wq to quantita- 
tively characterize the quality of crystalline structures. The local order pa- 
rameter that measures the degree of 6-fold-orientational ordering is defined 
gS] as 

n(k) 

Mn) = W) j:^ (3) 

where dkj is the angle between the vector r^ — Yj and a fixed direction in 
the plane (chosen here to be the positive x direction) and the sum over j 
extends over all n(k) nearest-neighbor of the fc-th atom as provided by the 
DT. The global order parameter associated with bond-orientational order is 
then obtained as an average over all particles 
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fe=i 

In a static ordered triangular solid, we have n(k) — 6 and 9k j is a multiple of 
7r/3 for all j = 1, ... , 6. In such a case, | (ifg) |= 1 and the plot of (i^g) on a 
complex plane will result in a single spot on a circumference with unit radius. 
Because of the fluctuations due to the zero-point motion we expect that, even 
in the perfect crystal, this spot will be broadened along the circumference 
and toward the center. In the presence of vacancies or other defects we expect 
significantly less orientational order and the plot of Wq should extend over a 
larger portion of the circle (in the limiting case of the liquid system, it turns 
in a distribution centered around the origin). 



3 Results: periodically replicated crystal 

Recently we have performed the first systematic study of multiple vacan- 
cies in two dimensional (2D) solid 4 He [35] . Dealing with low temperature 



50 



40 



* 30 - 



20 



10 - 



"" 


AE = 25.6 K + n v x 0.61 K 




ii 


AE = n v x 7 K 






- 


" /$! 






- 


/ s 




□ 


M = 240 - 


- /i 




A 


M = 480 - 


/ 




A 


M = 490 


j 




• 


M = 960 




. . . . I , . , , I , , , , 





M = 1440 " 



10 



15 



20 



25 



30 



35 



Fig. 1 Defect formation energy AE nv at constant density as a function of n v 
in 2D solid 4 He at p = 0.0765A -2 computed in boxes with different lattice site 
numbers M. Dashed lines are linear fit to the data. 



properties, 4 He atoms are described as structureless zero-spin Bosons inter- 
acting through the standard HFDHE2 Aziz potential [SO]. From the com- 
puted equation of state we found the freezing and the melting densities to 
be, respectively, p f = 0.0672 A" 2 and p m = 0.0724 A" 2 . In our SPIGS com- 
putation for the periodically repeated system we have used the pair-product 
approximation [51] for the imaginary time projector; as a compromise be- 
tween good accuracy and reasonable computational effort 5t has been set to 
1/40 K _1 . A total projection imaginary time r = 0.775 K" 1 ensures con- 
vergence to the ground state. We have simulated a 2D crystal at densities 
slightly above the melting density. The geometry of the SB is such that it 
is compatible with a regular triangular lattice and the study has been per- 
formed with M = 240, 480, 960 and 1440 lattice positions and PBC were 
applied in all directions. When N = M we have an ideal ordered crystal. By 
removing atoms from the initial ordered configuration we can study a system 
with a variable n v = M — N number of putative vacancies. n v ranged from 
to 30, with a maximum vacancy concentration x v — n v /N of 0.025. For 
larger values of x v the crystalline order is lost. After removing n v particles 
from the starting ideal ordered configuration the dimensions of the SB were 
rescaled to restore the original crystal density p. As a test of convergence we 
performed simulations starting from different positions of the n v vacancies. 

The crystalline order is stable for all the cases considered and the prop- 
erties of the system display a dependence on n v rather than on x v . Both 
the diagonal and the off-diagonal properties were found to exhibit a double 
regime behavior [35] . For n v < 6 the defects are easily identified in the snap- 
shots of the atomic positions and the defect formation energy AE nv is found 
to be monotonic and systematically sublinear with n v . This is consistent with 
the presence of an attractive interaction among vacancies (as confirmed by 
the vacancy- vacancy correlation function [35]). For n v > 10 AE nv becomes 
linear in n V) with a cost of 0.61K for each additional vacancy, i.e. the forma- 
tion energy for an additional vacancy in this regime costs about ten times 
less than the energy (7 K) required for a single vacancy (see Fig. [T]). A simi- 
lar double regime behavior has been found also in the static structure factor 



and in the amount of disorder Xd- For n v < 6 the main Bragg's peaks are 
found as sharp as in the perfect crystal but with a decreased height and Xd 
is found to increase with n v (as expected for an increasing disorder). When 
n v ^ 10 the behavior is quite different: the main Bragg's peaks are slightly 
broadened but with an integrated intensity that is almost independent of n v 
and also Xd does not depend on n v . By looking to the sampled configurations 
we realized that multiple vacancies for n v < 6 are highly correlated with a 
preference to form fluctuating linear structures. On the other hand when 
n v > 10 the vacancies inserted in the initial configuration lose their iden- 
tity and transform themselves into quantum dislocations [35 . In the range 
6-10 there is a crossover between the two regimes. We have now analyzed 
such multiple vacancy systems in terms of the orientational order parame- 
ter $6- $6 also displays a similar double regime behavior. Some examples of 
the obtained results are reported in Fig. [5] By comparing the result of ($%) 
in the ideal perfect crystal (Fig. [5^) with the results for the crystal with 1 
vacancy (Fig. [5Jd), 4 vacancies (Fig. [3fc) and 10 putative vacancies (which 
equilibrate into dislocations, Fig. [2jl), one can see that, in the cases when 
disorder is present, (&&) fills more of the complex plane far from the point 
5ft[(<Jg)]=l on the real axis with respect to the perfect crystal case. However, 
in going from the 4 vacancies to the dislocations case, the signal of disorder 
decreases, inverting the trend from 1 to few vacancies; again this is a sig- 
nature of the double regime behavior found in monitoring other quantities. 
An interstitial (Fig. [5^) gives a distinct signature of disorder with respect 
to that of a vacancy; 4 interstitials give a disorder map in (t?g) that is less 
extended than that for a single interstitial, showing some cooperative corre- 
lations also among interstitials. In fact, by looking at the DT of the sampled 
configurations we can recognize the presence of dislocations also in this case. 

This double regime is found also in the large distance behavior of the one- 
body density matrix p\. For small n v , as n v increases, the plateau in the pi 
tail, which is the signature of ODLRO, decreases, and this can be interpreted 
as an effect of the vacancy-vacancy interaction. Should this trend continue 
into the large n v regime, no ODLRO could be present in the thermodynamic 
limit, implying that an extended crystal with vacancies would not be super- 
solid. However we find a different behavior: for n v > 10 the plateau in the 
large distance tail of pi is restored (see Fig. [3]) mainly due to the ability of the 
dislocation cores to transfer particles among them. In fact, dislocations turn 
out to be very mobile and are able to induce exchanges of particles across the 
whole system, which is a necessary condition for ODLRO. Unfortunately, no 
size scaling analysis on p\ has been done yet in the dislocation regime because 
of the prohibitive computational cost of off-diagonal simulations of systems 
large enough to accommodate a large number of dislocations. In Fig. [3^ and 
Fig. [3}d one sees the presence of ridges in the tail of p\. As discussed in 
detail in Ref. [35] , these ridges are due to commensuration effects between 
dislocations and SB size. 

An important question is the behavior of many vacancies in a 3D crystal. 
Based on finite temperature simulations, statements in the literature [25] 
are found that in presence of many vacancies the crystal becomes unstable 
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Fig. 2 (Color online) Intensity plots of the probability of W'g on the complex plane 
for a bulk 2D 4 He crystal at p — 0.079 A -2 , for the (a) ideal perfect crystal 
with TV — 572, (b) crystal with a vacancy (JV = 571), (c) crystal with 4 vacancies 
(TV = 568) (d) crystal with dislocations (n„ = 10) with TV = 562, (e) crystal with an 
interstitial (TV = 573), (f) crystal with 4 interstitials (dislocations) with TV = 576. 
Note the logarithmic intensity scale. 





Fig. 3 One-body density matrix pi(r, r ') in the x — x' , y — y' plane computed in 
a defected 4 He 2D crystal (n„ = 10) at p — 0.0765A -2 for different numbers of 
lattice sites M: (a)M = 480 (b) M = 700 (c) M = 960 . 



against separation into a vacancy-rich and a crystalline vacancy-free phase. 
We found no such instability at T — K not only in 2D but also in 3D for 
as many as 98 vacancies as discussed in Ref. [35 . 



4 Periodic boundary conditions, the crystalline state and the 
issue of zero point defects in solid 4 He 



All simulations deal with a finite number N of particles and, in order that 
such finite size systems mimic a bulk system, PBC turn out to be very useful. 
In a fluid system the extrapolation to the bulk limit of the simulation results 
for short range properties is rather straightforward, unless one is close to a 
critical point, because the results depend smoothly on N. The situation is 
more complex for a crystalline solid because simulated properties are affected 
also by commensuration effects between the SB and the crystalline unit cell. 
As an example of this consider the following. The stable phase of solid 4 He at 
low T and pressure is the hep crystalline state. However, from exact quantum 
simulations of realistic models of 4 He it turns out that all three phases hep, 
fee and bec are stable (with a different value of the energy) for a suitable 
choice of the SB geometry and of N [23]; i.e. the average density, p, has 
to be large enough and the sides of the SB have to be integer multiples of 
the sides of the unit cell of the crystal under study. Finite size scaling for a 
given phase is performed by considering a special set of numbers of particles, 
for instance 180, 448, 900 for hep, but this explores a special set of states 
predetermined by the choice of SB and N. The conclusion is that finite size 
scaling is very special in a crystalline solid and it does not give an unbiased 
estimate of the properties of the bulk system. One might think that a way 
out of such a limitation is to start from a disordered configuration and let the 
system crystallize into the preferred state. The problem is that PBC again 
constrain the system. In fact, compatibility between a crystal lattice and 
PBC for a given SB is present only for certain orientations of the crystal axis 
with respect to the axis of the SB. The spontaneous crystallization of the 
particles has essentially zero probability to pick up such special orientations 
and, in fact, computer experiments |52| show that crystallization does indeed 
take place, but the crystal axis take arbitrary orientations with respect to 
the SB so that the crystalline order is deformed in order to comply with the 
PBC. 

For similar reasons one cannot expect to see during the computation the 
spontaneous appearance of a vacancy. An objection [25 , 12 , 26 to the presence 
of ground state vacancies is that such defects cost large energies, quantum 
simulations [18,21,25,25] agree on a value of order of 15 K in 3D at the 
melting density, so that no defect should be present at low T. This, however, 
fails to recognize that those simulations addressed specifically the question 
of the energy cost of a vacancy as an excited state and this is a separate 
question from the presence of vacancies in the ground state. It is instructive to 
consider the lattice gas model introduced by Mullin [27J . This model displays 
a liquid phase, a normal solid and a supersolid phase. In the supersolid state 
the occupation probability of a lattice site is less than unity and delocalized 
vacancies are present (more precisely vacancy-interstitial pairs) as an effect 
of zero point motion. Yet, the excitation spectrum has a gaped vacancy- 
interstitial branch. Anderson |15] has argued that in a Bose quantum crystal 
a unit population of a lattice site is incompatible with quantum mechanics 
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and that some vacancies have to be present in a highly quantum system like 
solid 4 He. 

One way to address the question of the presence of vacancies in the ground 
state of solid 4 He is to proceed in a way similar to what is done in classical sta- 
tistical mechanics in computing the equilibrium concentration of vacancies in 
a classical solid. The method is based on considerations of a macroscopic sys- 
tem, not the one that is simulated, and a basic assumption is that vacancies 
can be treated as weakly interacting objects. Within variational theory, using 
a SWF as ground state, a finite concentration x v of zero-point vacancies has 
been found with i„ ~ 1.4 x 10 -3 at melting density [53 . This is a sizeable 
concentration not too far from the actual experimental upper bound. Exten- 
sion of this computation to the case of the exact ground state is a difficult 
and open problem also because one has to take into account that vacancies 
are strongly interacting and tend to transform themselves into dislocations, 
at least in 2D, as discussed in the previous section. 

Another approach to study the presence of ground state vacancies is to 
get away from PBC by considering a confined system. This is the topic of 
the next section. 



5 Confined crystal 

Specific simulation details: In the study of confined systems we have opted 
for an implementation of the worm algorithm |llj with a fixed number of 
particles and for the pair-Suzuki approximation |44j for the imaginary time 

propagator e~ br¥L '; St has been set to 1/360 K" 1 and the total imaginary 
projection time r = 0.775 K _1 is equal to the value used in the simulations 
of the periodically replicated crystal. 

The confining potential has been chosen to be circular symmetric so that 
all directions are equivalent and, in addition, this ensures minimal mismatch 
with the crystalline triangular structure so that the interfacial region should 
be minimal. Since we are not interested in simulating any specific confine- 
ment, we have chosen the confining potential V cx t to be only repulsive and 
to represent an hypothetical surrounding infinite 4 He bulk crystal at density 
Poxt with the atoms considered as a continuum. Then V cx t reads 



T/ , , _ f V L j{R-r)-V L .i(R-r min ) R-r min <r<R 
V ^ r >-\0 r<R-r min 



(5) 



where R is the radius of the confining potential and Vlj is obtained by 
integrating a standard 12-6 Lennard- Jones potential on the half-plane: 



37T£CT 6 
VXj(r) = /9ext 



21cr 6 \ / 1 



80r l 



10 I J 



(6) 



and r m j n is the position of the minimum in Vlj . As e and a we use the typical 
Lennard- Jones parameters for 4 He, i.e. 10.22 K and 2.556 A respectively, 
and p C xt is set to 0.0765A -2 . Such a confining potential fails in accounting 
for the confining curvature, but the curvature effect is small for the large 
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Fig. 4 (solid line) Radial density profile, p(r), for a system with N — 187 and 
R — 40.0 A. (dotted line) Confining potential, V ex t{r) 



values of the radius R here considered. The initial configuration is obtained 
by building a regular triangular lattice at a density ps and discarding all 
the particles falling out the disk of radius R — r m j n - The number of particles 
inside this disk is referred to as iVo(ps). If ps is below a certain value the 
atoms quickly lose the initial crystalline order and the system evolves into a 
liquid state. As an example, in Fig. 0] the case of a confining potential with 
R = 40 A and N (p s ) = 187 (ps = 0.0435 A" 2 ) is shown. From the plot 
of the radial density profile it is evident that the confinement gives rise to 
an interfacial region that covers an annulus whose thickness is not larger 
than about 8 A. By increasing the number of particles at fixed radius R, we 
find that the particles become more and more localized and beyond a certain 
value of N the solid order becomes visually evident. For the solid phase, 
such interfacial region is expected to be thicker because of the effect of the 
mismatch between the triangular geometry of the confined crystal and the 
circular one of the confining potential. We have considered different values 
of the radius R and we find that the interfacial region never penetrates more 
that 15 A inside the system. 

In order to promote the presence of defects in the inner region, additional 
simulations are performed by subtracting or inserting particles at random 
positions from an equilibrated configuration for the system with Nq{ps) par- 
ticles. For the confined crystal the quantity Xd, introduced in the previous 
section, is no longer appropriate because of the lack of a conservation law 
like the one in ([T]). In place of Xd, as an estimate of the local disorder, we use 
the density pa of the particles that are miscoordinated, i.e. with coordination 
number different from 6. 

Results: We have considered three different values for R, 44.6 A, 54.6 
A and 64.6 A, that give rise, respectively, to Nq(ps) = 433, 685 and 931 4 He 
atoms when the starting configuration is obtained by cutting from an ideal 
triangular lattice at ps = 0.0765 A -2 . When equilibrated, the system is solid 
with an average density p av slightly different from ps because of relaxations 
of the atoms in the confining potential (/o av are reported in the legend of 
Fig. [6]) . In the present paper we discuss mainly the case of the intermediate 
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Fig. 5 (color online). Delaunay triangulation of typical configurations of a 2D 
confined 4 He crystal with N = 671 and 7? = 54.6 A (a) after few hundreds 
Monte Carlo steps from the random removal of 14 particles from an equilibrated 
configuration of the system with N = No(ps) and (b) when equilibrium has been 
reached. We report the coordination number for atoms only when different from 6. 
The couple 5-7 indicates a dislocation core in 2D and the clumps of two 5 and two 
7 indicate bound pairs of dislocations. 



radius value R = 54.6 A. In addition to the "magic" number N (ps) — 685 we 
have considered also an initially "defected" crystal by changing the number 
of particles. From an equilibrated configuration with 685 4 He atoms we have 
removed 14 particles at random positions obtaining a system with N — 671 
particles, and we have also randomly inserted 14 particles in order to obtain 
a system with N — 699 particles. In the case of such initial defected state 
the MC evolution has two stages. At first, after a few hundred MC steps, the 
crystal relaxes around the injected defects and we have a confined solid with 
defects (vacancies or interstitials) that are present also in the inner region. 
In Fig. [5] (a) the DT of one typical configuration in this initial evolution 
for the N = 671 system is shown. At a later stage of the MC evolution, 
the system reaches an equilibrium state with no sign of defects recorded in 
the inner region (see Fig. [S] (b)). This is confirmed also by the density of 
miscoordinated particles p*. From the plot in Fig. |5]it is easy to recognize in 
the system with R = 54.6 and 64.6 A a disk of about 30 A radius in which 
pa is essentially compatible with the expected value in periodically repeated 
ideal crystals at similar densities. In Figj6]yO0(r) for a periodically repeated 
crystal with one vacancy is also shown. Its value is about 40 times larger of 
that of the system without a vacancy. It is clear that the values of pd (r) in 
the inner region of the confined system are incompatible with the presence 
of such defect. 

The suggestion of a defect-free inner region is supported also by the fol- 
lowing argument. Consider a periodically repeated crystal. As discussed in 
Sec. [21 even in the ideal perfect crystal there is a finite probability to find 
miscoordinated particles due to the zero-point motion. Since this is a fluctu- 
ation effect, we might expect the permanence time ta, i.e. the mean number 
of MCS in which a given ill coordinated particle remains ill coordinated, is 
smaller in an ideal perfect crystal compared with the case of a crystal with a 
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Fig. 6 Local densities, pg{r), of particles with coordination different from 6 as 
function of radial distance for a number of 2D confined 4 He crystals. For comparison 
with the values of bulk 2D 4 He crystals at similar densities are also shown (see 
legend) . 



vacancy, because the vacancy is a permanent defect whose dynamics is gov- 
erned by the rate of vacancy jump. Indeed, at density p — 0.77 A~ 2 we find 
£0 = 4.0 MCS for the periodically repeated ideal perfect crystal and t$ = 9.5 
when a vacancy is present. A similar computation for the confined system 
and considering only particles with r < 30 A gives t$ = 3.0, 2.7 and 2.8 in 
the systems with N = 671, 685 and 699, respectively. Thus not only is the 
inner region free of defects, but it turns out to be somehow less fluctuating 
than the system with PBC, presumably due to the stabilizing effect of the 
confining potential. 

Also the study of the orientational order parameter ipg corroborates the 
finding that no defects permeate into the inner region of the confined crystal. 
From the plots Fig. [7] it is evident that within 30 A from the center the 
results for (t^) in the confined system with R = 54.6 are close to the results 
obtained in the periodically repeated ideal crystal reported in Fig. [2] for the 
three numbers of particles, N = 671, 685 and 699. In the range 30-40 A there 
is a crossover to the outer region of width ~ 15 A strongly perturbed by the 
confining potential. 

In conclusion, our study of a confined system shows that the inner region 
is free of defects like vacancies or interstitials; if defects are present in the 
initial configuration these anneal away and go to the interfacial region. Since 
our system is finite, our computation allows to give only an upper bound to 
the concentration of ground state defects. Taking into account that the inner 
region has a radius of 35 A in which the average number of 4 He atoms is 
about 300, this upper bound is 3 x 10~ 3 . We notice that this upper bound is 
above the estimated concentration of vacancies given by a SWF (in 3D) 53 . 



14 







n n n 



-0.5 0.5 

»[<#e)] 



-1 -0.5 0.5 



»[0&5>] 



-0.5 0.5 



Fig. 7 (Color online) Intensity plots of the probability of !^6 on the complex plane 
computed for different annular shells {r m in < r < r max ) for 2D confined 4 He 
crystals with R = 54.6 A. (column a) N = 671 (column b) N = N (p s ) = 685 
(column c) TV = 699; (first row) < r < 20 A (second row) 20 < r < 30 A (third 
row) 30 < r < 40 A (fourth row) 40 < r < 50 A (fifth row) 50 < r < R A. Note 
the logarithmic intensity scale. 
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6 Conclusions 

The issue of ground state defects in crystalline 4 He is still debated. In 3D 
we have studied up to 98 vacancies in system of 2548 atoms and in 2D up 
to 30 vacancies in a system of 1410 atoms and in all cases the crystalline 
order, as studied by different methods, turns out to be stable as long as the 
concentration of vacancies is below about 2.5%. Multiple vacancies are not 
weakly interacting though, they display signature of strong interaction. The 
most detailed study has been performed in 2D. For up to 6 vacancies one 
find a tendency for the vacancies to form linear structures. Starting from 10 
vacancies, in the relaxed state such vacancies cannot be identified anymore 
but the system displays the presence of fluctuating dislocations. This suggests 
that should ground state defects be present in solid 4 He, such defects would 
be dislocations and not vacancies. 

We have stressed the role of periodic boundary conditions in quantum 
simulations of a crystalline solid. The computed quantities are affected not 
only by the usual size effects, but also by commensuration effects between the 
simulation box and the crystalline lattice. We have investigated solid 4 He in 
2D without the use of PBC by confining the atoms via an external potential. 
We find that the effects of the confining potential extend mainly in a layer 
of width of order of 15 A, inside this interfacial region the atoms display a 
crystalline order similar to that of a bulk crystal and we find no evidence for 
the presence of defects. If defects like vacancies or interstitials are introduced 
in the initial state, such defects migrate to the interfacial region leaving the 
inner region free of defects. Taking into account the finite size of the system, 
this allows us to put an upper limit of 3 x 10~ 3 to the concentration of 
zero-point defects. 

Many extensions of the present study can be foreseen. A stronger bound 
on the concentration of zero-point defects will require the study of larger 
systems. The walls of a real container are never smooth and it will be inter- 
esting to see how the properties of the systems are affected by assuming a 
confining potential not as smooth as the one used in the present study. The 
study of confined system in 3D should be interesting but it is hindered by 
the very large number of particles needed in order to get relevant results. 
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